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Abstract 

This paper investigates a computational strategy for studying the interactions 
between multiple through-the-width delaminations and global or local buckling in 
composite laminates taking into account possible contact between the delaminated 
surfaces. In order to achieve an accurate prediction of the quasi-static response, 
a very refined discretization of the structure is required, leading to the resolution 
of very large and highly nonlinear numerical problems. In this paper, a nonlinear 
finite element formulation along with a parallel iterative scheme based on a multi- 
scale domain decomposition are used for the computation of 3D mesoscale mod- 
els. Previous works by the authors already dealt with the simulation of multiscale 
delamination assuming small perturbations. This paper presents the formulation 
used to include geometric nonlinearities into this existing multiscale framework 
and discusses the adaptations that need to be made to the iterative process in order 
to ensure the rapid convergence and the scalability of the method in the presence 
of buckling and delamination. These various adaptations are illustrated by simula- 
tions involving large numbers of DOFs. 

keywords: nonlinear multiscale computation; domain decomposition method; de- 
lamination; buckling; composites 

1 Introduction 

Delamination is one of the main degradation mechanisms of laminated composite ma- 
terials. This phenomenon is generally initiated by large interlaminar stresses due to 
edge effects, impacts, concentrated loads or macroscopic defects. Under some loading 
and geometric configurations (e.g. compression and a high slenderness coefficient), 
buckling is likely to occur once the delaminated zone has reached a critical extent dur- 
ing the propagation phase of the delamination process. Then, this geometric instability 
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can lead to an increase in interlaminar stresses, an acceleration of the delamination rate 
and, eventually, to the failure of the structure. The first analytical studies of buckling 
and delamination growth in the 70-80's were done by 1 15 8[|5][T3); in the last decade, 
new analytical studies were proposed by |[T2l|28l[T9ll . In a finite element context, the 
first works were based on fracture mechanics l30l [32l l27l while in more recent publi- 
cations cohesive models were also used to deal with geometrically nonlinear problems 
Il7irn i29l . Asymptotic numerical methods ITUl [171 l6l were also applied to delamina- 
tion buckling problems. Despite these many contributions to a better understanding of 
the mechanics of laminated composites, the inclination by industry to substitute virtual 
simulations for expensive experimental tests raises new issues. Thus, the numerical 
prediction of combined buckling and delamination remains a scientific challenge be- 
cause, even when using calculations on the mesoscale 11221 . a highly refined discretiza- 
tion of each ply is necessary in order to describe the delamination fronts and buckling 
loads properly. Therefore, multiscale and parallel computational techniques are being 
developed for buckling ifTTl |251 and debonding problems lfl4l [181 . 

In this work, we propose a mixed and multiscale domain decomposition strategy 
for the parallel simulation, in a geometrically nonlinear context, of composite lami- 
nates which are subject to multiple delaminations. Our approach was adapted to the 
treatment of geometric nonlinearities from an existing LATIN (LArge Time INcre- 
ment) multiscale strategy for delamination problems under the assumption of small 
perturbations lfT8l . 

Here, the geometrically nonlinear evolution is handled through a total Lagrangian 
formulation and delamination is modeled on the mesoscale using a cohesive interface 
model based on damage mechanics JT|. For this first-time approach, the intralaminar 
degradations are considered to be negligible and the layers are assumed to follow a 
hyperelastic law. Unilateral contact conditions are introduced by means of an interface 
law in order to avoid interpenetration over the delaminated surfaces. The reference 
problem and the substructuring process are summarized in Section[2] 

The LATIN strategy [20 1 consists in dividing the structure into volume substruc- 
tures separated by 2D interfaces, both of which are mechanical entities. As a result, the 
reference problem associated with the chosen mesomodel is naturally substructured, 
and both the unilateral contact conditions and the cohesive interfaces are handled at 
the interfaces of the domain decomposition. Section[3]introduces the LATIN algorithm 
proposed for the resolution of the nonlinear substructured problem. For the sake of 
computational efficiency, three scales are considered in this resolution: 

• The microscale corresponds to small-wavelength phenomena, which occur be- 
tween neighboring substructures. 

• The macroscale corresponds to the permanent verification of a weak form of 
equilibrium throughout the structure. This part of the algorithm, which makes the 
method scalable, is achieved through the definition of a small number of macro- 
scopic degrees of freedom per interface, which must satisfy continuity conditions 
and are linked together by a homogenized behavior constructed automatically. 

• In some cases, the number of substructures and interfaces (which depends on the 
number of plies) may be so large that the macroscopic problem cannot be ad- 
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dressed by direct solvers. Therefore, the substructures are grouped into "super- 
substructures" (whose size is determined by the available processor memory), 
and the macroscopic problem is solved using a primal domain decomposition 
method (25] . The third scale (or supermacroscopic problem) is introduced clas- 
sically in the course of balancing the supersubstructures with respect to the rigid 
body modes. 

This framework is very much under the control of the operator. In order to pilot the 
calculation of slender structures, the following parameters must be adjusted: 



The influence of neighboring subdomains and interfaces is represented by what 
is known as "search directions" (also called "interface impedances" in ET\ ). 
These parameters must be adapted to the aspect ratios of the slender structures 
by introducing well-chosen anisotropic coefficients; this point is discussed in 
Section |4~T1 



• Unilateral contact with or without friction between small surfaces can be handled 
successfully by the multiscale LATIN method ||9]|23l. However, as shown in Sec- 
tion 4.2 in the case of contact between slender structures over large delaminated 
areas, the search directions should be optimized according to the interface's state 
(open or closed) because incorrect values could generate artificial stiffnesses or 
induce interpenetration of the contact surfaces. 

• Because of the stiffness loss as a result of buckling and delamination, the macros- 
tiffness and search directions may become irrelevant and the macroscopic opera- 
tors may need to be adjusted in order for the homogenized behavior to represent 
the current state of the structure better, as illustrated in Sections [4. 1.2| and |4.3| 

• The supermacroscopic problem is solved using a projected preconditioned conju- 
gate gradient algorithm, which requires the setting up of a convergence threshold. 
The solutions to this problem were developed in [18] and will not be repeated 
here. 

With these improvements, the multiscale analysis of large combined buckling and 
delamination problems becomes possible. The capabilities of the strategy are illus- 
trated by two examples involving geometric instabilities coupled with debonding (Sec- 
tion!^]). 



2 The reference problem 
2.1 Notations and assumptions 

Figure [T] shows the general motion of a deformable body. The body is considered to 
be an assembly of material particles M identified by their initial coordinates X_ with 
respect to the Cartesian frame Bo = {X 1 ,X 2 ,X 3 }. In general, the current positions 
of these particles are defined by their coordinates x with respect to another Cartesian 
frame B — {x^x^^x^}. ^ n tn i s work, the two coordinate systems Bo and B are the 
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Figure 1 : The general motion of a deformable body 



same, but we will refer to them as separate entities in order to associate each quantity 
with the initial or the current configuration. The displacement u of particle M between 
the two configurations is defined as: 

u = x-X. (1) 
The deformation gradient tensor F is given by: 

£ = I x = I d + £ M, (2) 

where V Q denotes the gradient with respect to the initial configuration. Let us note 
that F, which transforms vectors in the initial configuration into vectors in the current 
configuration, is called a two-point tensor. 

Let us consider a laminated composite structure E (see Figure[2} occupying domain 
O bounded by dVl in the current configuration C, and consisting of Np adjacent plies 
P. Each ply P, of mass density pp, occupies a domain Sip such that £7 = {J PeE fip. 
The plies are assumed to be separated by Np — 1 cohesive interfaces. The structure 
is subjected to an external surface traction field F d over part dVlp d of the boundary 
<9£1 and to a displacement field U_ d ) over the complementary part dfljj d - The body 
force per unit mass is denoted f d . A ply P defined in domain Sip is connected to an 
adjacent ply P' through an interface Tpp> = dilp n dflp'. Let T = {J PeE Tp, where 
Tp = Up'eE Tpp'- The relevant quantities of E (e.g. volume, area, surface tractions 
or density) can be described in reference to the configuration before deformation Co- 
The index o will be used to denote the initial (undeformed) configuration, e.g. flo, 
dflo, T p P „F d p Po . 




Figure 2: The reference problem: undeformed and deformed configurations 
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The objective of the present work is to study the response of E subjected to a pre- 
scribed loading starting from the initial configuration, and resulting in large displace- 
ments and rotations accompanied by progressive damage of the interfaces T. Because 
of the very small thickness of each layer of the composite (« 0.125 mm), the delami- 
nated areas can be very slender {L de i ammated /h ply > 100). 

Our study of this problem relies on the following assumptions: 

1. structure E may undergo large displacements; 

2. the behavior of the plies is hyperelastic; 

3. the loads are independent of the configuration of E {i.e. follower forces are not 
considered); 

4. the evolution over time is considered to be quasi-static; 

5. isothermal conditions are assumed; 

6. the displacements along T can be large, but the displacement discontinuities in 
the non-fully delaminated part of the structure or in the contact region are small; 

7. r is assigned irreversible softening behavior by means of an interface law con- 
necting tractions with displacement discontinuities; 

8. in the delaminated region, corresponding points in adjacent plies may separate, 
regain contact after separation, or remain in contact. 

Assumption|6]enables the displacement of interface T pp, to be defined as the mean 
value of the displacements of the adjacent plies (see Figure [3J: 

(u PP ,) = i(u P , + up) , over r PP / . (3) 
Thus, interface Tpp' itself is defined as: 

r PP > = {x\x = X + (u PP ,) ; over T PP ,} . (4) 

With this definition, Tp p, is viewed as the mean surface between deformed plies P and 
P', which, in fact, due to Assumption [6j can be considered to coincide geometrically 
(i.e. same area and same orientation). Because of this geometric coincidence, the 
interface's deformation gradient tensor can be defined as: 

Z PP , - Y (u PP >) + L d , over Tpp, . (5) 

Thus, interface Tp p> can be viewed as a zero-thickness medium that carries out the 
transfer of traction forces between plies. The displacement gap of interface Tpp, is 
given by: 

[up P /] = Up, — Up , over Tpp, . (6) 
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2.2 The nonlinear damage interface law 



An anisotropic interface law can be formulated as a relation between the traction vector 
t P and the displacement discontinuity vector [uppi] : 

t P = KppXhpp'}) kp P <] , over r PP , , (7) 

in which a dependence on the current orientation of the interface Tp pi is introduced. 
One can show, as proven in [1 1, that law (|7]l satisfies the principle of material frame 
indifference if As sumption [6] is satisfied. 

The expression of the local stiffness operator E_ ppi of interface Tpp> can be made 
explicit with respect to a local orthonormal frame B n = {n 1 ,n 2 ,n 3 } moving together 
with the interface, where n 3 is the deformed unit normal to T ppi pointing from P to 
P' (see Figure |5J: 



Kpp^K^n^n^ (8) 



where: 



(1 - di) k° t 

(Kl° cal )=[ (l-da)*? 

(l-h + ([u P p,}-n 3 )d 3 )k° n 



h+ denotes the positive indicator function. fc° and are the initial elastic stiffnesses 
of the interface, with the dimension of a force per volume. The softening behavior 
of the interface model when the structure is loaded is simulated by the introduction 
of the dimensionless scalar damage variables with values ranging from (healthy 
interface point) to 1 (completely damaged interface point). We use the evolution law 
defined in |]4), which has the advantage of using a single damage variable to handle 
different macroscopic delamination modes of the interface. 



Pi 



X lA Po 




Figure 3: The normals to the undeformed and deformed interfaces 

To calculate the internal power corresponding to the interfaces, it is mandatory to 
express the traction vector as a function of the displacement discontinuity vector with 
respect to frame Bo- In order to do that, one needs to write K pp , = Kf^X^ ® X y 
Introducing the orthogonal matrix which characterizes the transition from B n to Bo'- 

Q = Q l3 X. l ®n J Vi,j = {1,2,3}, (9) 
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where Qij — 2Li 1 Rj are the direction cosines of X_ t relative to n,-, the interface's local 
stiffness becomes: 

(K?* obal ) = Q ip K l p ° cal Q jq Vi,j = {l,2,3}Vp ) « = {l,2,3} (10) 

(with implied summation over indexes p, g). 

From Nanson's formula [j] which describes how an infinitesimal surface element 
deforms during a given motion, the unit normal vector at each point of the deformed 
interface Tpp> is calculated as follows: 



(11) 



-PP' — 3 



where 7V 3 is the initial (undeformed) unit normal to T p p , pointing from P to Pg (see 

■0*0 

Figure [3j. 

Finally, one needs to calculate the traction vector in the initial configuration: 

tPo = Jr^P = det{g ppl ) \\g- pl Nz\\ tp . (12) 

After an interface has been fully damaged, a frictionless contact law is assumed: 
the gap between two plies must remain nonnegative, and only compression can be 
transmitted when the plies are in contact. 

Remark: The relative interface displacements are assumed to remain small enough 
(Assumption |6]l for the contact to be detected only between points which were con- 
nected in the initial configuration. 

2.3 Substructured formulation 




Figure 4: Substructuring of the laminated composite structure 

The laminated composite structure E is divided into substructures and interfaces 
as shown in Figure [4] Each of these mechanical entities possesses its own kinematic 
and static unknown fields related by its constitutive law. The substructuring process is 
governed by the objective of matching the domain decomposition interfaces with the 
cohesive material interfaces, so that each substructure belongs to a unique ply P and 
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its behavior is geometrically nonlinear. In the initial configuration Co, a substructure 
E defined in domain Q Eg is connected to an adjacent substructure E' through an 
interface r„ „, = dft e H dH, E > (see Figure |5). The surface entity r„ „, applies 

force distributions F_ E , F E i and displacement distributions W_ E , WLe' to an< ^ E'o 
respectively. Let T Eo = Ue^e T e o e> - 

For a substructure E such that T Eo n (dClp d U dCljj d ) ^ 0, the boundary condition 
(Fd , U do ) is applied through a boundary interface T Ed ■ 




Figure 5: Unknown fields on the interfaces and substructures 
Let u„ be the displacement field, E „ the Green-Lagrange strain tensor, E „ the 

—E r —E ° ° —E 

Lagrangian strain rate, n „ the second Piola-Kirchhoff stress tensor, F _ the deforma- 
tion gradient tensor and J E — det(F ) the Jacobian of the motion at each point of 
substructure £7 . At each point of interface r„ _,, the displacement field is defined as 

Eo-Eq 

W P , the displacement gap as [W„ ] = W E , — W_„ , the mean value as (W„ ) 
h(W p, + W „ ) and the deformation gradient tensor as F „„, = VAW„ ) + I ,. 

£ x E Q ' ° —EE' =0 x E ' —d 

Then, the substructured quasi-static problem consists, at each step of the time in- 
tegration, in finding s — (s„ )„ „, where s p = (u p ,n ,W_ P ,F„ ) verifies the 
following equations: 

• Mass conservation of substructure E : 



Pe Je = PE , over Vt Eo (13) 

• Nonlinear kinematic admissibility of substructure E : 

Ee = \ faPE, + %^E + %^E ^E ) > ° Ver ^Eo d4) 

^|an^^|r E / 0Verr ^ (15) 

• Global nonlinear equilibrium of substructure E : 

V(u Eo , W Eo ) eU° E xW E , such that u% ]d n Eo = 

I k Eo : I«) dn = I PEg f_ d ■ u Eo dn + I F Eg - W% dV , 

(16) 
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where E(u%) = ^u Eo + %u Eo + %u E ^u Eo + 

• Hyperelastic orthotropic behavior of substructure E : 

=Eo = dE~ ' OV6r Eo ' ( ?) 

— Eo 

where ip is the stored energy function or elastic potential per unit of undeformed 
volume. In this first study, we use ib = \\L E , E_„ '■ E„ ■ 

z — Eo — Eo 

• Constitutive equation of interface T „ „, : 

E E 

1l E E ,([W E ],F E ,F E , F ) = , over Y E E , e T Eo . (18) 

• Behavior of the interface at the boundary V E do ' 

K Ed0 (W Eq ,F E J = , over 1^ . (19) 

The formal relation 1Z„ „. — 0, called the "interface behavior" and defined over 

E E - 

r„ can be made explicit in the three cases we are concerned with: 

E E 



• Perfect interface: 



• Cohesive interface: 



App'(\W E ] , F , F EEf 



(20) 



(21) 



• Unilateral contact interface (without friction): 

F F + F E , = 

— E„ — ^o — 



n 3 ■ [W E ] > and n 3 ■ F E >0 
(n 3 -[W Eg ])(n s -F Eo )=0 ° 
PF p = PF E , = ° 

— E„ — a n — 



(22) 



where substructures E and E' belong respectively to plies P and P' ; operator 
App' was introduced in Section 2.2 n 3 was defined in Eq. ([TTJ and P is the 
corresponding tangential projection operator. 

As already mentioned at the end of Section [2~2] contact in the delaminated inter- 
faces is detected only between points which had the same initial position. 

Remark: Usually, the admissibility of forces is written in the deformed config- 
uration; here, it can be easily converted to the initial configuration thanks to the 
geometric coincidence of the deformed interface (Assumption^. 
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Figure 6: The linear macrobasis for a plane interface 



3 The numerical resolution strategy 

3.1 The macroscopic scale 

In order to ensure the scalability of the method, one can solve a coarse global linear 
problem. The definition of the macroscopic fields required to formulate this problem 
refers only to the interface's unknowns. The action/reaction principle F E + F^ E i = 
is verified regardless of the interface's behavior. The objective of the macroscopic 
problem is to ensure that part of this equation is verified at any time: 

/ (E Eo + F E ,) ■ W M * dT = 0, VW M * £ W M , (23) 

where the displacement macrospace W M and its dual space F M are parameters of 
the method. These subspaces are common to neighboring substructures and induce a 
separation of the interface quantities which is made unique by the uncoupling of the 
virtual works: 

V(F Eo> W E J G T Eo x We, F_ Eg = F^+F™ , W Eq = W™ + , 

/ F Eo ■ W Eo dr = J • W% dr + J • W% dV , (24) 

Tfi ^e 

Usually, one chooses a common basis for the kinematic and static macroscopic 
fields of the interface. Numerical tests have shown that in order to ensure the numerical 
scalability of the method the macroscopic basis should extract at least the linear part 
of the interface forces (see Figure [§}. Indeed, this macroscopic space contains the part 
of the interface fields with the longest wavelength. Consequently, according to Saint- 
Venant's principle, the micro complement has only local influence. 

3.2 The iterative algorithm 

In this section, the iterative LATIN algorithm, which enables one to deal with nonlinear 
problems, is adapted to the resolution of the geometrically nonlinear substructured 
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Figure 7: Schematic representation of the LATIN iterative algorithm 

reference problem with nonlinearities localized at the interfaces. The finite element 
method is used to discretize the equations. 

The equations of the problem are divided into two groups: 

• Admissibility of the substructures and macroscopic admissibility of the inter- 
faces: 

- mass conservation of the substructures, Eq. ( fT3l >; 

- nonlinear kinematic admissibility of the substructures, Eq. ( p"4}[T5] i; 

- nonlinear static admissibility of the substructures, Eq. ( p"6| ); 

- behavior of the substructures, Eq. ( fT7| ; 

- macroscopic admissibility of the interfaces (after the linearization of the 
previous equations), Eq. ( |23j ). 

• Local (nonlinear) equations at the interfaces: 



interface behavior, Eq. (18 19 1 



The interface solutions s = (s p )„ „ — (W„ ,F„ ) „ cT? of the first set of equa- 
te £ '0 fcllj ^0 £/ fclli ^ 

tions belong to space A<i, while the interface solutions J = Cs„ )_ _ = (W_„ , F W )„ _ 
of the second set of equations belong to T. The converged interface solution s re f is 
such that: 

s ref € A d n r . (25) 

The resolution process consists in seeking the interface solution s re f alternatively in 
these two spaces: first, one finds a solution s n in Ad, then a solution s~ n+ i in F. In 
order for the two problems to be well-posed, one introduces the search directions E + 
and E~ which link the solutions s and s during the iterative process (see Figure |7b. 
Thus, an iteration of the LATIN algorithm consists of two nonlinear stages which are 
described in detail below. 

Remark: The LATIN approach [20, 24 1 is a general computational strategy for the 
resolution of time-dependent nonlinear problems which operates over the entire space- 
time domain. In our case, time is irrelevant and the capabilities of the LATIN method 
are not fully exploited. The LATIN strategy is also based on the idea of separating the 
difficulties and dealing with global linear equations and local nonlinear equations inde- 
pendently. In this work, we consider a global nonlinear stage (called the admissibility 
stage) which is solved as a series of global linear equations. 
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3.2.1 The local stage: 

In the local stage, the following local problems are solved at each point of interfaces 

r 

E E 

f n EoE ,([w E },F Eg7 F E ,,l EoE ,) = o 

Find (F E M Eo ,l E ' ,W E ' )™ch that: J (l Eo - F Eg ) - fc+ - W Eg ) = 

(26) 

The last two equations of this system define the search direction E + . In the case of a 
cohesive interface, Problem ( p6*| l is nonlinear and is solved using a modified Newton- 
Raphson algorithm. 

The following local linear problems are also solved at each point of boundary in- 
terfaces r E do ■ 

- — f K-E d (W w ,F„ ) = 

H„d ( F E ,W E J such *c ■ _ ■ _• fi (27) 

^ n » iS fi » K » 



3.2.2 The admissibility stage: 

The admissibility stage consists in solving the nonlinear equations in each substruc- 
ture under the constraint of macroscopic admissibility of the interfaces. This stage is 
carried out using an iterative Newton-Raphson procedure; therefore, at each iteration, 
the macroscopic admissibility is prescribed on the linearized system of equations in the 
substructures. 

In absence of macroscopic admissibility, the search direction E~ which couples 
the interface displacement and the force fields from the linear stage is: 

(F E -F ) + k~ (W -W E )=0, overly, . (28) 
The monoscale version of the nonlinear problem in each substructure becomes: 

V(u Eo ,W Eo )eU E xW E , 

I K Eo E (u ) : E(u Eo ) dn + I k- W . W Eo dT 
Jn Eo a ° Jt Eo e o 
s v ' 

-P lnt {u ,w ) 

= I Pe L ■ u Eo dn+ f (F Eg + k E W ) ■ W Eo dV . (29) 

S v ' 

This problem is solved using a Newton-Raphson algorithm. At each iteration i, the 
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correction term (J'Su^ , % 5W ^ ) is calculated by linearizing d29): 

Jn Eg 

+ [ K £o g{ l 8u E ) : dtt + f k~ *6W E -W Eo dT = P ext +P mt ( l u E , *W 



(30) 



where EM = H^+'V^+'^'m^ ¥jL+ t ¥j)Vjuw ), ) is known 

from the previous iteration, and , t+1 W_ w ) = ^Vl w ) + ( 4< 5m,t /<^W" f ). 

Then, at each iteration, macroscopic admissibility is applied to linearized problem 
pOj l. The satisfaction of the macroscopic equilibrium of interface forces ( |23j ) suffices 
to ensure the scalability of the method. 

Condition ( |23] l is incompatible with search direction ( |28] i; hence this search di- 
rection is weakened and verified "as best can be" under the macroscopic constraint. 
Technically, this is achieved by using a Lagrangian multiplier whose stationarity leads 
to a modified local search direction: 



VW% € W° , [ ( t+1 F E - F ) • W% dVo 

/■pi 



k- Eo C +1 w Eo -w E )-k E ^w yw Eo dr = o. (3i) 

In order to solve this system, a relation linking l + 1 F 1 ^ and 2+1 W_„ is derived 

from the subdomain equilibrium plus the modified search direction ( |30||3l| l. The linear 
problem to be solved in order to find ( t 6u p , l SW w ) becomes: 

^0 ^0 

y(u*E ,W* Eo )eU E xW E) 

I K Eo E E fu Eo ) : (%*Su Eo Y u Eo +%u Eo Y;6u Eo ) dQ 

+ [ K Eo K{ l 5u E ) : E(u* Eo ) dn + [ k- l 6W E -W% dT = l P ext +P in t( i u E , i W E ), 
Jn Eo Jt Eo a <> oo 

(32) 



f ■^i~ M 

where % P ext = P ex t + / k l+1 W_ ■ W Eg dT . One can prove that if K Eo and 

Jr Eo E o 

k are symmetric, positive definite operators, then Eq. ([32) is well-defined and has 

E o 

a unique solution. Due to the linearity of Eq. (|32), one can define a linear relation 
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between the interface displacements and the loading: 

VZ^eTl, / i+1 W E ■ F Eg dT 

i Eo (l Eo + k E J +1 W%) + ^W Eo +^W E ) ■ F* Eg dT , (33) 



where F E = F_ p + k W_ p , and operator l H„ is the dual Schur complement of 

E o * o 

substructure E modified by the search direction at iteration i (which depends on the 

geometric configuration of the previous iteration i), while % W C results from the con- 



densation of the volume loading and of Pi n t( l u„ , l WL^ ) onto interface IV.. l W ' 
are the interface displacements of the iteration i. 

The corresponding interface forces are obtained using the modified search direction 
( f3T| and projected onto the macroscopic space: 

i+i F . w M * dr = ( (<L« l+1 W_ M + 1+1 F E ) • W M * dT , (34) 

where: 

VW M * e W M , 



f %™ l+1 w M ■ w M *dv = [ (k~ - k- k- y +1 w M ■ W M *dT 

Jr Eo E o J Teo e e e e 

' X Ke -W M *dr = [ (l Eo ~k E (*M E l Eo VW E + i W E ))-W M *dT . 

It.. 



Ve 



l h is classically viewed as the homogenized behavior of substructure E and is cal- 

E o 

culated explicitly for each substructure by solving local subproblems ( |32| ), taking the 
vectors of the macroscopic basis as boundary conditions over T E . Clearly, l h M de- 

E o 

pends on the geometric configuration of the previous iteration. 

Finally, Relation ( |3"4"| > is introduced into Eq. ( |23~| , which expresses the admissibility 
of the macroforces, leading to what is called the macroscopic problem: 

~M* x f . ...~M ~M* 

f ~M* t-^ f ■ , , ~M ~M* 

= E / ^o-W dTo-Y, + Fe -W dT (35) 
Jdn F , Jr E „ 

The macroscopic problem is discrete by nature. Therefore, it has an algebraic form 

l L M l+l W_ = i+1 F M , where l+1 W_ is the vector of the components of the La- 
grange multiplier in the macroscopic basis. 
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The right-hand side of Eq. ( |35| l can be interpreted as a macroscopic static resid- 
ual obtained from the calculation of a monoscale linear stage. In order to derive this 

~M 

term, Problem ( |32] i must be solved independently in each substructure, setting l+L W 
to zero. The resolution of the macroscopic problem of Eq. ( |35j ) leads to the global 

■4-1 ~ M 

knowledge of Lagrange multiplier l+L W_ , which is finally used as a prescribed dis- 
placement to solve the substructure-independent problems of Eq. ( [32] ). 

The resolution of ( [32] > in the substructures is carried out using the finite element 
method. 



Remarks: 



• The convergence of the algorithm can be improved by introducing a relaxation 
stage after the admissibility stage. The admissibility solution s n is renamed s„; 
then the relaxed solution s n is defined by: 

s„ = ^s„ + (1 - M)«n-i j (36) 
where /i is a relaxation parameter usually taken equal to 0.8. 

• The LATIN error indicator adopted here was successfully used in Ifl8l l2ll. This 
criterion is based on a measure of the non-satisfaction of the constitutive laws of 



Eq. (18 19 1 in the admissibility stage, since these are the only equations which 
are not verified at this stage. More precisely, each time a solution s„ £ Ad 
is obtained, an indicator of the convergence of the algorithm is calculated by 
integrating the local residuals of the interface behavior over the structure. 



4 Analysis of the parameters of the algorithm in the 
case of slender structures 

The search direction parameters of the local stage (Af^ gE ^ and of the admissibility 
stage (k~ )( E are symmetric, positive definite operators which represent the influ- 
ence of the neighboring subdomains and interfaces, such as interface impedances or 
Schur complements of the rest of the structure fio \ CIe - It was empirically shown in 
previous studies [ 9 , 2 1 1 that there is an optimum set of these operators which depends 
on the interface behavior. 

For the monoscale approach (which does not include the satisfaction of the macro- 
scopic problem) applied to massive isotropic and homogenous structures with contact 
or perfect interfaces, the basic setting is the scalar approximation k + = k~ ~ E/Lq, 

where E is Young's modulus and Lq is a characteristic length of the structure. The 
multiscale approach enables most of the effect of the search direction to be localized, 
so the classical setting becomes k + = k~ ~ E I Lr„ , where Lr„ is a characteristic 

to E o E 1 E o E o 

length of the interface. 

In the case of cohesive interfaces with damage, the local search directions k + of 

the cohesive interfaces must be set to infinite values, because values which are too 
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small can lead to the stagnation or the divergence of the algorithm, as explained in 
lfl8l . This choice also enables the interface's quantities to be calculated directly in the 
local stage, so local Newton iterations are unnecessary. Concerning search direction 
k~ , the optimum value would be twice the actual interface stiffness 2k°(l — d) (k° 

E o 

denoting the undamaged interface's local stiffness), which would be equivalent to pre- 
scribing the exact interface behavior as an interface condition in the admissibility stage 
fl8l . Unfortunately, the use of this value would require the operators to be updated 
very frequently, which would be expensive; instead, a monitoring strategy has been 
proposed in lfl8l . 

Our particular study, which focuses on multiple buckling-delamination interactions 
in composite laminates, is characterized by the occurrence of bending in very slender 
geometries because the plies are very thin (« 0.125 mm) and small delaminated ar- 
eas can have high slenderness coefficients (L delaminated I 'h p i y ^> 100). As a result, 
classical values for massive structures are completely inadequate and the following 
difficulties need to be tackled: 

• Concerning the search direction for perfect interfaces, the scalar E/Ly Eo nearly 
equals the stiffness of the neighbors (disregarding the geometry of the structure), 
which could affect performance in the case of slender structures, especially in 



bending and in buckling (see Section 4. 1 . 1 



• The occurrence of buckling leads to major changes in the deformed configuration 
which require additional continuity conditions in order to ensure convergence, as 



proven in Section 4. 1 .2 



• In the case of multiple buckling-delamination interactions, it is possible for de- 
laminated surfaces to separate, regain contact or remain in contact throughout the 
evolution, so during the calculations the same surface can go through different 
states which are unpredictable. Unfortunately, in the case of slender structures 
and subdomains, there may be different suitable search directions depending on 
the state (open or closed), and an inadequate value could lead to nonphysical so- 



lutions or to stagnation of the iterative process. See Section 4.2 for the proposed 
remedies. 

The use of optimum search directions k~ for cohesive interfaces lfl8l requires 

the operators to be updated and reconstructed according to the damage state at 
each point of the cohesive interface, which leads to an increase in CPU time. 
In the procedure proposed in 1 18], the search directions are updated only when 
an interface becomes fully damaged. Unfortunately, the efficiency of this proce- 
dure deteriorates when the size of the substructures increases, and becomes even 
worse in the case of slender subdomains with large cohesive interfaces (because 
the updating is performed only once the entire interface is damaged, which can 



take several time steps). See Section 4.3 for the recommended solution 



Finding the most favorable values requires not only additional identifications in 
order to take into account the geometry of the structure, but also the introduction of a 
new step: the micro/macro separation of the search direction. 
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If one divides search direction E of the admissibility stage into a macroscopic 
part E~ M and a microscopic part E~ m , Eq. ( pi) can be rewritten as: 

Vf M *6W M , 

(F- -F E ). W M * dTo + [ k~ M (W - W_ -W_l). W M * dT = 



T-i — u -°0 /t^ E, 

1 En •> I 



(37) 



J T E E J r E E E 



where W m is the space orthogonal to W with respect to the inner product L (Te ). 
Parameter k~ M represents the stiffness of the interface for the macroscopic prob- 

lem. Its optimum value, if it exists, is the homogenized interface behavior: for perfect 
interfaces (infinite stiffness), this parameter must be as large as possible; for homoge- 
nous elastic interfaces {e.g. cohesive interfaces with constant damage), it is related to 
the current stiffness k~ M = 2fc°(l — d) (where k° is the undamaged interface's local 

E o 

stiffness). Parameter k~ m is the micro part of the search direction and, classically, 

E o 

k~ m = E/Lr E (k~ rn = 2k°(l — d)) is considered to be a good starting value for 

perfect interfaces (elastic interfaces). 

In order to deal with the difficulties one at a time, we will seek the optimum pa- 
rameters independently for the different types of interface behavior using academic 
examples of 3D slender structures. In the following sections, separate analyses will be 
presented for perfect interfaces in bending and buckling (Section 4. 1 1, contact inter- 
faces (Section pk2] ) and cohesive interfaces (Section |4~3] >. 

All the analyses were carried out using a fully parallel C++ program, taking ad- 
vantage of the three scales of the mixed domain decomposition method proposed. The 
transfers of data among the different processors required for the parallel computations 
were performed using the MPI libraries. Each processor was assigned to a set of con- 
nected substructures and their interfaces and used to calculate the associated operators 
and solve the local problems. This was achieved technically thanks to a METIS routine 
and helped reduce the number of interfaces among the processors. 

4.1 Perfect interfaces 

Here, we present a study of the search directions for slender structures. The analysis 
concerns a cantilever plate in bending under the assumption of small perturbations 
and shows how a thin geometry affects the convergence rate and the scalability of the 



strategy (Section 4.1.1 1. The proposed improvements as a result of the bending study 



will be tested in the case of a buckling example in Section 4. 1 .2 
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4.1.1 Bending: 

The problem being considered is that of a slender plate which is built-in along one side 
and subjected to a bending load _F d in the form of a surface force distribution along the 
opposite side, as shown in Figure [8] The material is isotropic and homogeneous, the 
geometry was fixed to: Lq = 20 mm, ho = 0.2 mm and bo — 1 mm. The substruc- 
turation was made only along the X\ -direction, modifying the number of substructures 
to study the convergence rate. Figure [8] shows the substructuration in 8 subdomains. 
The whole mesh totaled 1.8 million DOFs with 10 linear wedge elements through the 




Figure 8: The 3D cantilever plate in bending: geometry, substructures and interfaces 
Table [T]presents the number of iterations as a function of the number of substruc- 



tures. Row A corresponds to the use of the classic values k 



-M 



E/Lt 



while Row B corresponds to continuous macrodisplacements k — s- oo, keeping 

E o 

k~ m = E/Ly e ■ In these two cases the convergence rate depends on the number of 

substructures and deteriorates when this number increases, even if the continuity of the 
macrodisplacements is enforced. 





number of substructures 


g 


16 


32 


64 




number of LATIN iterations until convergence (0.1%) 










A 


without continuity of the macrodisplacements 
and without anisotropic search directions 


15 


15 


22 


38 




number of LATIN iterations until convergence (0.1%) 










B 


with continuity of the macrodisplacements 
and without anisotropic search directions 


14 


15 


21 


36 




number of LATIN iterations until convergence (0.1%) 










C 


without continuity of the macrodisplacements 
and with anisotropic search directions 


3 


6 


9 


15 




number of LATIN iterations until convergence (0.1%) 










D 


with continuity of the macrodisplacements 
and with anisotropic search directions 


4 


4 


4 


5 



Table 1: Influence of the number of substructures on the convergence rate (Lo/ho = 
100) 

We propose an enhancement to the microscopic directions based on the fact that 
the local stresses and displacements in the normal and transverse directions are very 
different. Plate theory leads to the following relations between the orders of magnitude 
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of the normal and shear stresses and between the orders of magnitude of the normal 
and transverse displacements: 



O 



T x 1 x 3 
h 



O 



(J x 1 x 1 
L 



O 



h 



O 



ux 3 
L 



Therefore, since the search directions relate tractions to displacements, they should 
take anisotropic values: 



(39) 



where the ratio Lo/ho is a macroscopic quantity representing the whole structure, the 
normal value of (k~ m ) n — E/Lr E remaining constant. 



Using these improved search directions and setting k = k ' 

E„ E„ 



the number of 



iterations dropped by more than 60% (Row C of Table[T]i, and the strategy became even 
more efficient and scalable when the continuity of the macrodisplacements k~ AI — >• oo 

was enforced (Row D of Tableyj. 
4.1.2 Buckling: 

In this section, we consider the geometrically nonlinear example of a plate built-in at 
both ends and subjected progressively to a negative end displacement, with a perturba- 
tion consisting of a central force. The data are: Lo — 10 mm, ho = 0.1 mm, b = 1 
mm, E = 135,000 MPa and v — 0.3. The geometry was divided into 640 substruc- 
tures and 1,464 perfect interfaces (see Figure |9j. The mesh totaled 2.2 million DOFs 
with 12 linear wedge elements through the thickness. The macroscopic problem and 
the supermacroscopic problem represented 13,176 DOFs and 372 DOFs respectively 
and were solved using 64 processors. 

The nonlinear buckling analysis was performed in 96 time steps, resulting in the 
evolution of the axial compression load as a function of the transverse displacement 
of the central perturbed point shown in Figure 10 This numerical response agrees 



perfectly with the theoretical response given in B3TII . which has a critical Euler force 




undeformed 
configuration 



Figure 9: The initial configuration and the final deformed shape after the last time step 
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4.4 N. Figure [9] also shows the final deformed configuration. 



equal to P c = 

It is important to note that buckling began close to the 6 th time step 

A search direction study similar to that presented before for bending under small 
perturbations was carried out for this buckling problem. Figure 1 1 shows the num- 
bers of LATIN iterations for each of the first 20 time steps corresponding to three 
different settings. In the first case, corresponding to the application of the values 

k~ M = k~ m = E/Lr E , convergence was not reached after the 3 rd time step. Using 

e e •- 

the same values, but this time with anisotropic differentiation of the search direction, 
the method diverged before buckling (after the 6 th time step). Convergence could be 
achieved only by increasing the value k~ M = k~ m = E/Lr E at least a hundredfold 

(which is close to ensuring the continuity of the macrodisplacements, see the A-bars 
in Figure 111. The continuity of the macrodisplacements (k" 



-M 



00) seemed to be a 



necessary condition for convergence in the case of large displacements (see the B-bars 
in Figure 111. With this macroscopic continuity, the anisotropic differentiation of the 



search directions improved the convergence rate (see the C-bars in Figure 1 1 



Remark: The satisfaction of the macroscopic problem in the admissibility stage (Eq. 
( [35) l) plays a major role in the handling of global geometric nonlinearities thanks to the 
transmission of the large-wavelength part of the solution. This requires the expensive 
continuous updating and assembling of the macroscopic homogenized operator l L M , 
which depends on the current configuration of the substructures. There have been some 
unsuccessful attempts to update the macroscopic problem only at the beginning of each 
time step, or even only after each local stage, but they resulted in divergence problems 
or erroneous solutions depending on the time step discretizations. In our work, we used 
continuous updating, but were able to limit the updating and, thus, the computation 
time thanks to the introduction of some additional criteria. 



4.2 Contact with possible opening 

Usually, in contact problems, no micro/macro separation of the search directions of 
the admissibility stage is applied (k~ AI = k~ m = k~ ) and the value fc + = k~ — 

Eg E E E E 

E/Lr Eg seems to be a good choice [23 1. In the case of open interfaces in slender struc- 
tures, the best values are k + = k~ = 0. The mechanical behavior of the structure is 

so different with and without contact that an incorrect setting of the search directions 
can affect the convergence of the method dramatically. Therefore, we had to develop 
a search direction updating strategy which takes into account the predicted status of 
the interface. This updating strategy must be a compromise between performance and 
stability. 

Then we applied this updating strategy to an opening contact problem and a closing 
contact problem, both initialized with an incorrectly predicted interface status, in order 
to verify whether the updating algorithm leads to the correct solution. 
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4.2.1 An opening contact interface: 

The first example concerns a plate built-in at both ends with an initial central delamina- 
tion ao, subjected progressively to a negative end displacement U_ d and to a perturbation 



consisting of a symmetrical central force F_ d , as illustrated in Figure 12 The data are: 



Lq = 20 mm, h = 0.2 mm, bg = 1 mm and ao = 10 mm. The geometry was divided 



into 80 subdomains and 126 interfaces (see Figure 12 1 



Figure 13 shows the evolution of the error in the first time step as a function of the 



number of LATIN iterations. Using the correct guess fc + = k~ = 0, convergence was 

E o E o 

achieved in 10 iterations. Conversely, using incorrect values k + — k~ — E/Lr E , 

it was extremely difficult to obtain the correct solution: because of the additional in- 
terface stiffness, stagnation in a non-physical configuration occurred. Therefore, we 
undertook to check the physical status of the interface every 10 iterations and to update 



the search directions accordingly. Figure 14 shows different states obtained thanks to 
this updating strategy. 

4.2.2 A closing contact interface: 

The second example concerns a plate built-in at both ends with an initial central de- 
lamination ag subjected to a central vertical bending force F_ d applied to the lower ply, 



as illustrated in Figure 15 The data are: Lq = 20 mm, ho = 0.2 mm, bo = 1 mm and 



ao = 16 mm. The geometry was divided into 80 substructures and 126 interfaces (see 



Figure 15 i 



In that configuration, the optimum guess for initializing search direction k = 

E„ 



16 1. The introduc- 



k =E/Ly e led to a reasonable number of iterations (see Figure 

E o 

tion of an incorrect stiffness (assuming an open interface k + = k~ = 0) resulted in 

E o E o 

an erroneous configuration with penetration. It was necessary to update the status in 
order to obtain a suitable value of the search direction and evaluate the error criterion 



correctly. Figure 17 shows the deformed configuration before and after updating. One 
can observe that overlapping of the plies occurred at the beginning, but was eliminated 
after updating. 

Remark: The use of search directions k + = k~ = E/Lr E for the closed inter- 

E o E o 

faces and fc + = k~ =0 for the separated interfaces resulted in a proper macroscopic 

E o E o 

problem (representing the stiffness of the contact interface) and a correct converged 
solution. The proposed updating scheme according to the interface's state appeared to 
work properly for an interface going from a closed state to an open state or vice versa. 
However, in the case of more complex problems, the strategy was not always found to 
converge because it is difficult to find the exact setup of the updating algorithm in order 
to avoid divergence. After some empirical tests, we decided to use for both states a uni- 
fied search direction equal to k + = k~ = (E/Lr E )/(Lo/ho) 2 , where Lo/ho is the 

E o E o 

slenderness coefficient of a ply. This choice leads to physically sound solutions with 
only a few more iterations than using the optimum value for each state of the interface. 
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4.3 Cohesive interfaces 



Cohesive interfaces lead to the same difficulties as contact interfaces: the recommended 
values are k + — > oo and k~ M = k~ m = k~ = 2fc°(l — d), and the main problem 

E E E E 

resides in the definition of a stable and efficient updating strategy for k~ with respect 

E 

to the evolution of d. 

The strategies were evaluated based on a DCB test under small perturbations. Fig- 
ure 18 shows the deformed specimen with an original crack ao subjected to a vertical 
displacement U_ d separating the two arms formed by the crack and built-in at the op- 
posite end. The data are: Lo = 20 mm, ho = 0.5 mm, bo = 2 mm, ao = 10 mm, 
E = 135,000 MPa, v = 0.3 and k° = 100,000 N/mm 3 , a = 1, n = 0.5 and 
Y c = 0.4 N/mm. The structure was divided into 160 substructures and 324 interfaces 
(see Figure 18 1, which amounted to 550,000 DOFs for the whole mesh and 2,916 DOFs 
for the macroscopic problem. Because of its small size, the macroscopic problem was 
solved using a direct solver. There were at least 10 elements in the process zone and 
the load was applied in 16 time steps. 

Figure 19 shows the load-displacement curve at the end of one of the two arms of 
the specimen. The global response until the complete failure of the specimen agrees 
with the theoretical solution calculated for a beam on an elastic foundation |16| using 
linear fracture mechanics theory for the propagation 1 3 1 . The difference in slope which 
can be observed before the softening part of the curve is due to the fact that the dam- 
age law used in the calculation enables damage to take place progressively before the 
rupture of a point. For the comparison of the different updating procedures, only the 
first 16 time steps were considered. 

Figure 20 shows the crack front after the 15 th and 16 th time steps. It is important to 
note that a first element becomes fully damaged after the 10 th time step; subsequently, 
the process zone expands to a new subdomain with each new time step to reach five 
completely damaged interfaces after the 16 th step. 

Figure 21 shows the convergence rate of each time step and the corresponding total 
CPU time for the different updating schemes. As expected, the number of iterations 
increases after the 10 th step. Without updating, the number of iterations explodes (A- 
bars), although the CPU time remains moderate. If, as proposed in 11811 . updating is 
performed only at the interfaces which have been fully damaged (B-bars), both the 
number of iterations and the CPU time decrease respectively by a factor of 3 and a 
factor of 2. The optimum strategy seems to consist in updating at each Gauss point after 
a predetermined number of LATIN iterations (e.g. 100), in which case both the number 
of iterations and the CPU time decrease another 65% compared to the previous strategy 
and become very small. The more aggressive strategy which consists in updating at 
each Gauss point after each iteration (D-bars) leads to the smallest number of iterations, 
but a huge CPU time (25 times that of the procedure C). 



22 



5 Examples of combined buckling and delamination 



5.1 A built-in beam under compressive loading 

The example discussed here consists in the simulation of a built-in plate with a central 
initial delamination ao subjected to an axial compressive loading U d and a symmetrical 



central perturbation F_ d (see Figure 22 1. The data are: Lo = 20 mm, ho = 0.2 mm, 
b = 1 mm, a = 10 mm, E = 135,000 MPa, v = 0.3, k° n = k° t = 100,000 
N/mm 3 , a = 1, n = 0.5 and Y c = 0.4 N/mm. The geometry was divided into 



1,280 substructures and 3,248 interfaces (see Figure 22 1, leading to a mesh totaling 
1.4 million DOFs with 24 linear wedge elements through the thickness (12 elements 
in each ply). The macroscopic problem and the supermacroscopic problem contained 
respectively 29,232 DOFs and 132 DOFs, and were solved using 24 processors. 

The solution of this problem is characterized by the competition between global 
and local buckling. The critical local buckling load is approximately twice that of a 
beam of length ao built-in at both ends: 



2 ( 47r ^ ;/ °) = 7IU()N ' 



where Io is the second moment of area of a single ply. The global buckling load of a 
beam with an initial crack is given approximately by the formula JT]: 

pfobal = (Lo^o) f^ESIo) , ( 2ao\ (^EJo\ = 



Lo J V J \ Lo J \ Lq 

In the case of local buckling without imperfection, the propagation condition can 
be approximated by the formula Q: 

3 tf4 , oe2 _ 4G c b 



16 ? ^ ^plocal ' 

where £ = WL / 2 /a>o is the dimensionless transverse displacement parameter and G c 
is the critical energy release rate (for the damage interface law used here, G c = Y c ). 
Thus, the propagation condition in local buckling is: 

wlo/2 = 0.34 mm . (40) 



First, we studied the numerical response of the structure under a small symmetrical 
perturbation (F_ d = 0.2 N). From the evolution of the compressive load as a function 
of the maximum transverse displacement of the upper ply shown in Figure [23] one can 
see that before the load reaches the critical local value p l c ocal = 71.06 N the response 
switches from local buckling to global buckling (after the 14*' 1 time step), which re- 
duces the load to about the critical global value pa lobal = 44.41 N and induces mode- 
II delamination at about 0.75 mm maximum transverse displacement. The deformed 
configuration and the corresponding crack front after the last time step are shown in 
Figure 24 An increase in the amplitude of the symmetrical perturbation to F_ d = 2 N 
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is enough to determine the local buckling mode. Figure 23 indicates a local buckling 
load equal to about 56 N, which is less than the approximate evaluation. Mode-I de- 
lamination begins at about 0.36 mm maximum transverse displacement, which is very 
close to the theoretical delamination condition d40b. The deformed configuration and 



its crack front after the last time step are shown in Figure 22 The solution obtained can 
be compared to the numerical solution of a beam given in [ 1 1. There is good agreement 
in terms of critical load and delamination propagation. 

5.2 Multiple through-the-width delaminations 

This section deals with a built-in laminate made of 4 plies of thickness ho with two 
initial central delaminations oo at 2hg and 3ho, subjected progressively to a negative 
end displacement U_ d and a central perturbation F d applied to the upper ply as shown 



in Figure 25 The data are: Lq = 20 mm, ho = 0.2 mm, bo — 1 mm, ao = 10 mm, 
Ei = 185, 500 MPa, E 2 = E 3 = 9, 900 MPa, i/ 12 = i/ 13 = 0.34, j/ 23 = 0.5, 
G12 = Gi 3 = 6, 160 MPa, G 23 = 3,080 MPa, k° n = k° t = 100,000 N/mm 3 , a = 
1, n = 0.5 and Y c = 0.4 N/mm. The lay-up sequence is [0/90] s . The geometry 
was divided into 1,280 substructures and 3,064 interfaces (see Figure |25]>, leading to 
a mesh totaling 2 million DOFs with 12 linear wedge elements through the thickness 
of each ply. The macroscopic problem and the supermacroscopic problem contained 
respectively 27,576 DOFs and 168 DOFs, and were solved using 30 processors. 

A first computation without cohesive interfaces was carried out in order to calculate 



the post-buckling response of the laminate. In Figure 26 the solid lines represent the 
compressive load as a function of the transverse displacement in the middle of each of 
the three layers determined in the laminate by the two initial cracks. One can see that, 
first, the upper layer undergoes nonsymmetrical local buckling after 80 N. Then, at 100 
N, the middle and lower layers (i.e. the unperturbed layers) buckle in the opposite di- 
rection, leading to a symmetrical local buckled configuration. For the nonsymmetrical 
buckled shape, the transverse displacement of the middle layer is positive (pulled by 
the upper ply) and is separated from the lower layer (point A). Subsequently, in the 
symmetrically buckled shape, the transverse displacement of the middle ply becomes 
negative (pulled by the buckling of the lower layer) and the lower and middle layers 
are in contact (point B). Another calculation without initial delamination resulted in a 
global buckling load more than five times the critical load with cracks. 

The response taking into account delamination between plies is also shown in Fig- 



ure 26 The buckling behavior is almost the same as the response without cohesive 
interfaces. The critical buckling load is smaller because of the damage interface law. 
After point B, the load decreases due to the fact that the first elements in the crack 
front are completely damaged. At point C, the middle and lower layers are no longer 
in contact. The deformed configurations and crack fronts after time steps A, B and C 



of Figure 26 are shown in Figure 27 
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6 Conclusion 



We presented an efficient calculation of extensive delamination in the presence of ge- 
ometrically nonlinear effects thanks to a three-scale domain decomposition strategy 
based on an iterative algorithm. This method enables the resolution of huge nonlinear 
systems of equations on the most suitable scales. 

The LATIN algorithm we proposed can handle both the geometric nonlinearities 
in the substructures and the surface degradations at the interfaces, thus making the 
treatment of both nonlinearities easier. The second scale (i.e. the macroscopic problem) 
enables the rapid transmission of the large-wavelength part of the response through the 
introduction of an updated homogenized behavior of the substructures which takes into 
account the deformed configuration and the current damage state of the interfaces. The 
macroscopic problem was solved using a parallel iterative solver. 

The classical values of the method's parameters (i.e. the search directions) were 
shown to be inadequate for slender structures. Therefore, we modified them in order 
to ensure the scalability of the method, its independence with respect to the geometry 
of the subdomains, an efficient convergence rate and an adequate CPU time for the 
treatment of combined buckling and delamination. Examples showing the capabilities 
of the strategy were also presented. 

In subsequent developments, the introduction of large sliding contact conditions 
in the delaminated area and the coupling of this 3D model with a plate model in the 
low-gradient zones should be envisaged. 
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Figure 10: The load-displacement curve for a compressed cantilever plate 
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Figure 1 1 : Influence of the search directions in the case of buckling 




Figure 12: The separation of a contact interface 
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Figure 13: The LATIN error in the open contact (first time step) for different search 
directions 
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Figure 14: Deformed configurations after updating (magnification factor x500) 




Figure 15: The closing contact interface 
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Figure 16: The LATIN error in the closing contact (first time step) for different search 
directions 




Figure 17: Deformed configuration after updating 




Figure 18: The substructures and interfaces of the deformed DCB specimen 
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Figure 19: The load-displacement curve of the DCB test 
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Figure 20: The crack front and the interfaces after the 15 and 16 time steps 
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Figure 21: The convergence rate and the corresponding CPU time 



32 




Figure 22: A compressed plate with an initial delamination in local buckling: substruc- 
tures, deformed configuration (with no amplification) and crack front 
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Figure 23: The load-displacement curve of the compressed plate with an initial delam- 
ination 




Figure 24: The compressed plate with an initial delamination in global buckling: sub- 
structures, deformed solution in global buckling (without amplification) and crack front 
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Figure 26: The load-displacement curves of the three layers of a compressed laminate 
with two initial delaminations. Note: only some time steps are represented by markers. 




Figure 27: The deformed configurations and crack fronts after three time steps of a 
compressed laminate with two initial delaminations 
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